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Negative streamer ionization fronts in nitrogen under normal conditions are investigated both in a 
particle model and in a fluid model in local field approximation. The parameter functions for the fluid 
model are derived from swarm experiments in the particle model. The front structure on the inner 
scale is investigated in a ID setting, allowing reasonable run-time and memory consumption and high 
numerical accuracy without introducing super-particles. If the reduced electric field immediately 
before the front is < 50 kV/(cm bar), solutions of fluid and particle model agree very well. If the field 
increases up to 200 kV/(cm bar), the solutions of particle and fluid model deviate, in particular, the 
ionization level behind the front becomes up to 60 % higher in the particle model while the velocity 
is rather insensitive. Particle and fluid model deviate because electrons with high energies do not 
yet fully run away from the front, but are somewhat ahead. This leads to increasing ionization rates 
in the particle model at the very tip of the front. The energy overshoot of electrons in the leading 
edge of the front actually agrees quantitatively with the energy overshoot in the leading edge of an 
electron swarm or avalanche in the same electric field. 
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I. INTRODUCTION 



Streamers (T), S| are growing filaments of weakly ionized 
non-stationary plasma produced by a sharp ionization 
front that propagates into non-ionized matter. Streamers 
are used in industrial applications such as lighting [E Q 
or gas and water purification [E S[, and they occur in 
natural processes as well such as lightning 0, H, or 



transient luminous events in the upper atmosphere [1C 
Therefore accurate modeling and simulation of streamers 
is of high interest. 

Most streamer models (see e.g. 0EI1I3I1 
EE EE EE EE SE S3) are so called fluid models for 
the densities of different particle species in the discharge. 
These models build on the assumption of local equilib- 
rium: transport and reaction coefficients in the conti- 
nuity equations are functions of local parameters only. 
If this parameter is the local electric field, we refer to 
this assumption as the local field approximation. This 
assumption is commonly considered to be valid as long 
as equilibration length or time scales are much smaller 
than the spatial or temporal gradients in the electric field. 
For the strongly varying electric fields within a streamer 
ionization front, the validity of the local field approxima- 
tion was investigated in [IE HE S3 • The general 'senti- 
ment' in these studies is that the approximation suffices 
for practical purposes and that more detailed methods 
tracking the behavior of individual particles lead to just 
minor corrections. 

Another recent result supporting the fluid approxi- 
mation for streamers was that even streamer branch- 
ing [2E SE S3 can be understood in terms of an in- 
herent instability of the fully deterministic fluid equa- 
tions [IE SI HE HE & These studies have shown that 
a streamer in nitrogen can reach a state in which the 



width of the space charge layer that creates the field en- 
hancement at the streamer tip, is much smaller than the 
streamer diameters; the streamer then can branch spon- 
taneously due to a Laplacian interfacial instability. 

However, despite success and progress of fluid approx- 
imations and simulations for streamers, there are three 
major reasons to reinvestigate the local field approxima- 
tion: 

1. Not all streamers are alike. Experiments as well as 
simulations show that rapidly applied high electric 
voltages can create streamers that are more than an 
order of magnitude faster and wider than streamers 
at lower voltages 31] . Whether earlier findings on 
streamers in lower potentials apply to those fast 
and wide streamers as well has to be investigated. 

2. The detection of x-rays emanating from lightning 
strokes [H, SE S3, SE| indicates that electrons 
can gain very high energies within early stages of 
the lightning event. Therefore runaway electrons 
within streamer and leader processes should be in- 
vestigated. Runaway electrons by definition violate 
a local approximation. 

3. Streamer branching is an inherent instability of a 
fully deterministic fluid model. However, fluctua- 
tions of particle densities might trigger this insta- 
bility earlier than they would occur in the fully de- 
terministic fluid model. In particular, in the lead- 
ing edge of an ionization front, particle densities 
are very low and the fluid approximation eventually 
breaks down. As the front velocity of this so-called 
pulled front [IE EE Sfi| is determined precisely in 
the leading edge region, single particle dynamics 
and fluctuations should be accounted for. 
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These three observations motivate our present reinves- 
tigation of the local field approximation for streamers. 
The starting point is a Monte Carlo model for the motion 
of single free electrons in nitrogen. We note that com- 
plete streamers have been simulated with Monte Carlo 
particle models before [33], however, a drawback of such 
models is that the computation time grows with the num- 
ber of particles and eventually exceeds the CPU space of 
any computer. This difficulty is counteracted by using 
superparticles carrying charge and mass of many physi- 
cal particles; but superparticles in turn create unphysical 
fluctuations and stochastic heating. 

In the present paper, we compare the results of a Monte 
Carlo particle model and a fluid model. We circumvent 
the problems caused either by a too large particle number 
or by the introduction of superparticles by investigating 
a small, essentially one-dimensional section of the ioniza- 
tion front as illustrated in Fig. [TJ We suppress effects of 
lateral boundaries by periodic boundary conditions. As 
the electric field essentially does not deviate from the pla- 
nar geometry within the region where the particle den- 
sities vary rapidly, a planar ionization front [l7l . |38| is 
a very good approximation of this inner structure. Of 
course, a planar front will not incorporate the electric 
field enhancement caused by a curved front [HI HI] , but 
this outer scale problem concerns only the electric field 
and can be dealt with through an inner-outer matching 
procedure (3f| [4(1 El| . Planar fronts allow us to investi- 
gate individual particle kinetics and fluctuations within 
the front and its specific strong spatio-temporal gradi- 
ents in a systematic way and within reasonable comput- 
ing time. 

In this paper, we concentrate on negative streamer 
fronts in pure nitrogen under normal conditions. We 
thoroughly discuss the case where the reduced elec- 
tric field at the streamer tip is 100 kV/(cm bar), and 
we summarize results for fields ranging from 50 to 
200 kV/(cm bar). The paper is organized as follows. 
In section [Til nrs t our Monte Carlo particle code and 
its numerical implementation are described. Then the 
derivation of the fluid model is recalled, and the numer- 
ical implementation of the fluid model is summarized. 
Then swarm or avalanche experiments in a fixed field are 
performed in the particle model, the approach of elec- 
trons to a steady state velocity distribution is investi- 
gated, and the parameter functions for the fluid model 
are generated. This sets the stage for a quantitative com- 
parison of front solutions in particle and fluid model in 
section Mil Here first the setup of planar front simula- 
tions is described, then the results of the planar front 
simulations within fluid and particle model and analyt- 
ical results are presented and compared. The emphasis 
lies on front profile, front velocity and ionization level 
behind the front. It will be shown that differences can 
be attributed to the electron kinetics in the leading edge 
of the front where the electric field does not vary, and 
that the electron energy distribution there agrees quan- 
titatively with that in the leading edge of an ionization 




Figure 1: (Color online) The relation between the full 
streamer problem and the planar fronts described in this pa- 
per: the left picture shows the narrow space charge layer sur- 
rounding the streamer head [l9l. l29l. |30| ; the width of the layer 
is much smaller than the streamer diameter which creates the 
characteristic field enhancement ahead and field suppression 
behind the front. The right picture shows a zoom into the in- 
ner structure of the space charge layer with an essentially pla- 
nar ionization front as treated in this paper. In the transversal 
direction, periodic boundary conditions are applied. 

avalanche or swarm. Section [IVl contains our conclusions 
on the validity of the fluid approximation. An appendix 
contains analytical approximations on the ionization level 
behind an ionization front. 



II. SET-UP OF PARTICLE MODEL AND FLUID 
MODEL IN LOCAL FIELD APPROXIMATION 

In this section, we summarize features of particle and 
fluid models, their numerical implemention and mutual 
relation as a basis for the quantitative comparison of ion- 
ization fronts in particle and fluid model in section iLTTl 

Our starting point is a model that contains all mi- 
croscopic physical mechanisms that are thought to be 
relevant for the propagation of a negative impact ion- 
ization front in pure nitrogen. It models the generation 
and motion of single free electrons and positive ions in 
the neutral background gas. While propagating freely, 
the electrons follow a deterministic trajectory according 
to Newton's law. The collision of electrons with neutral 
molecules is treated as a stochastic Monte Carlo process. 
Because the mobility of the positive ions is two orders 
of magnitude smaller than that of the electrons, ions are 
treated as immobile within the short time scales inves- 
tigated in this paper. Neutral molecules are assumed to 
have a uniform density with a Maxwellian velocity dis- 
tribution. The electron-neutral collisions, including all 
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relevant elastic, excitation, and ionization collisions, are 
treated with the Monte Carlo method. Electron-electron 
or electron-ion processes as well as density changes of the 
neutral gas are neglected as the degree of ionization stays 
below 10~ 5 even at atmospheric pressure [H, HtJ. This 
well-known model will be summarized in section III Al 
The space charges can change the local electric field, 
this is accounted for by solving the Poisson equation. 
The particle model gives a very detailed and complete 
description at the expense of significant computational 
costs where we stress that one particle is one electron 
and superparticles are not used. 

If densities are high enough and fields vary slowly in 
space and time, the average behavior of the electrons can 
be modeled by a fluid approximation for electron and ion 
densities whose parameters depend on the local electric 
field only. The derivation of the fluid approximation can 
be formalized by taking the first and the second moment 
of the Boltzmann equation as is recalled in section Hi B 11 
However, for the practical purpose of determining mobil- 
ities, ionization rates and diffusion coefficients as a func- 
tion of the electric field, we directly perform swarm ex- 
periments with the particle model in a constant electric 
field. This procedure together with the averaging pro- 
cesses involved are described in section III CI Here also 
the relaxation of an electron swarm to steady state mo- 
tion and the velocity distribution of steady state motion 
in a given homogeneous field are discussed. 



A. The Monte Carlo particle model 

1. Physical processes 

In the particle scheme, at each instant of time i, there 
is a total number of N e (t) electrons and N p (t) ions. The 
single electrons are numbered by i = 1, . . . , N e (t) at time 
t; they are characterized by a position Xi(t) and a velocity 
Vi (t) , each within a continuous three dimensional vector 
space. Between collisions, electrons are accelerated and 
advanced according to the equation of motion 



m— =eE(x i ,i), — 



(1) 



where m and e are electron mass and charge and E(x, i) 
is the local electric field. For the positive ions, only their 
position x^ 1 , j = 1, . . . , N p (t) is taken into account while 
their mobility is so low that their velocity can be ne- 
glected. The electric field is determined from the Poisson 
equation 



V • E = 



e (n p 



(2) 
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Figure 2: Collision frequencies in N2 at 1 bar from [42 



respectively. Within a particle model, these densities are 



JV e (i) 



JV„(t) 

s (x,t) = ]T <5 3 (x-x^)), n p (x,i) = ]T <S 3 (x-x^)), 

i=i 

(3) 



where S is Dirac's 5 function. 



together with appropriate boundary conditions on the 
electrodes. Here n e and n p are electron and ion density, 



The collisions account for the impact of free elec- 
trons on neutral nitrogen molecules. As the neutrals are 
abundant, their density determines the probability of an 
electron-neutral collision. The collision can be elastic, 
inelastic or ionizing. In inelastic collisions, electron en- 
ergy is partially converted into molecular excitation; in 
ionizing collisions, electron energy is consumed to split 
the neutral into a positive ion and a second free elec- 
tron. The probability distribution of the different col- 
lision processes depends on the electron energy at the 
moment of impact; we use the cross section data from 
the Siglo Database [13] • Fig. [5] shows the frequencies of 
different electron-neutral collision processes for electrons 
as a function of energies in pure N2 at 1 bar pressure. As 
the collisions are random within a probability distribu- 
tion, the actual occurrence of a specific collision within a 
sample is determined by a Monte Carlo process. 

Once an elastic or inelastic collision process is chosen, 
the energy loss of the electron and therefore the absolute 
value of its velocity after the collision is fixed. However, 
model results will depend on the angular distribution of 
the emitted electrons, which again follows a probability 
distribution. Different scattering methods have been dis- 
cussed in the literature [H, [3, Hi, IH, E3| ■ Here we will 
only focus on the scattering method used in the present 
paper. 

In an elastic collision, the longitudinal scattering angle 
X and the azimuthal scattering angle <p relative to the 
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direction of the incident electron, are given as [44 

2pi(l-e(e)) 



cosx 



1 l + e(e)(l-2 Pl )- 



(4) 
(5) 



where the random numbers p± and P2 are drawn from 
a uniform distribution between and 1, and £(e) is a 
function of the incident particle energy e. The elastic 
cross section data [45[ for nitrogen was fitted empirically 
in Hi with 



0.065 e/eV + 0.26 Je/eV uJe/eV 

= ; ; (6) 

1 + 0.05 e/eV + 0.2^/e/eV 1 + AO^/eJeV 

In an inelastic collision, the energy loss of incident elec- 
trons has to be taken into account, but the scattering 
angle is calculated in the same way as for an elastic col- 
lision. 

In an ionizing collision, energy conservation dictates 



£l +£2 



(7) 



where e c , e% and £2 are the energy of the incident, the 
scattered and the ejected electron, respectively, and £i on 
is the ionization threshold energy. We use Opal's empir- 
ical fit for the distribution of the ejected electron energy 

H 



e 2 = B tan 



Pz arctan ■ 



2B 



(8) 



where B « 13 eV in the energy range of interest and p$ 
again is a random number equally distributed between 
and 1. For the scattering angles, Boeuf and Marode [46| 
assumed that (i) the incident, ejected and scattered elec- 
tron velocities are coplanar, and (ii) that the scattered 
and ejected electron velocities are perpendicular. These 
assumptions lead to 



cos xi = 



£1 



COS X2 



£2 



(9) 



where xi,2 are the respective scattering angles. The set 
of equations ([7]), ([5]), and © determines the distribution 
of energies and scattering angles of the scattered and the 
ejected electron in an ionizing collision. 



2. Numerical implementation 

The particle code moves electrons within the applied 
plus the self- induced field and includes their collisions. 
Therefore the numerical calculation consists of three 
parts: the Newtonian electron motion within the field, 
the field generated by the charged particles, and colli- 
sions. At each time step of length At, the field is calcu- 
lated from the charge densities on a lattice with grid size 
A£. Then the electrons move in continuous phase space 
according to the field, possibly interrupted by Monte 



Carlo collision processes. Electrons can experience more 
than one collision during one time step At. 

In more detail, position and velocity of the electrons 
are determined in continuous phase space from their 
Newtonian equation of motion ([!]) according to the elec- 
tric field at their initial position within the time inter- 
val. The commonly used integration is the leap-frog 
method (49| . in which the electron positions and veloci- 
ties are offset in time by At/2. 

For the electron-neutral collisions, time, type and scat- 
tering angles are determined in a Monte Carlo process 
by sequences of random numbers. For the data shown 
in Fig. [2 the maximal collision frequency z^ max is about 
9.7- 10 12 /s , therefore the minimal average collision time 
Tmin = 1/^max is about 0.1 ps. By introducing so-called 
null-collisions, in which no collisions occur, T m i n can be 
chosen as average collision time independently of the elec- 
tron energy e. The probability P(t) that an electron will 
travel a time t without collision (including null collisions) 
is 



P(t) = e-»» 



(10) 



Therefore the next collision time At co m S i on of an electron 
is drawn in a Monte Carlo process from the distribution 



At 



collision 



In- 



1 

P4 



(11) 



where p^ is again a random number. When a collision oc- 
curs, the energy of the incident electron is calculated, and 
the distribution of the collision processes is determined 
according to the collision frequencies from Fig. [21 then a 
random number determines the collision type (null col- 
lision, elastic, excitation or ionization collision). At the 
collision, the electron velocities are changed according to 
the processes discussed in Sec. Ill A 11 Then the next col- 
lision time for the particle is determined. This approach 
is described in more detail in [H, [H3| ■ 

At each time step of length At, the electric field is cal- 
culated on the grid with mesh A£. First, the number of 
elementary charges n p — n e within a grid cell is counted; 
it directly determines the charge density within the cell. 
Then the change of electric field components normal to 
the cell faces are determined from the densities within 
the cells through the Poisson equation This simple 
interpolation on cells of appropriate size causes no arti- 
facts as we are dealing with particles carrying just one 
elementary charge e, not with superparticles. The con- 
dition on the cell size is (i) that it is large enough that a 
single charge in the cell center does not create substan- 
tial fields on the cell boundary, and (ii) that it is small 
enough that no strong density gradients occur between 
neighboring cells. Here it should be noted that density 
gradients due to particle number fluctuations are strongly 
suppressed when we deal with real particles, not super- 
particles. Therefore more involved interpolation methods 
like Particle in Cell (PIC) are not required. 

The choice of the spatial and temporal mesh deter- 
mines the computational accuracy as well as the com- 
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putational costs. We have tested different meshes in 
planar fronts as described in section ITlI Al The results, 
most prominently the ionization density behind the front, 
converge for a sufficient discretization. However, a bal- 
ance has to be found between computational accuracy 
and computational costs. We choose the time step as 
At = 0.3 ps, which is of the same order as the min- 
imal average ionization time T m j n , and the cell size as 
A£ = 2.3 fxm, which is the basic length scale according 
to dimensional analysis in [HI]. On this mesh, the ion- 
ization density behind planar fronts has an error of less 
than 0.2%. 



B. The fluid model 

1. Derivation from the Boltzmann equation 

Fluid models in general are derived from the Boltz- 
mann equation [Ell, [H, For each species k, this 
equation describes the temporal development of the 6- 
dimensional phase space density /fc(x, v). For the elec- 
trons, omitting the subscript k, the Boltzmann equation 
reads 

| + Vl -v/-V„-f/=(f) ( . ( 12 , 

eE/m is the acceleration of the electrons under the influ- 
ence of the electric field. The term on the right hand side 
signifies the changes in the distribution function due to 
collisions that happen instantaneously. It is an integral 
expression that contains the phase space density /(x, v) 
both of the electrons and of the collision partners. 

Equation (TT2"|) can also be represented as an infinite se- 
ries of velocity moments. For the Z th moment, one multi- 
plies the Boltzmann equation (TT2"]) with an / th -order vec- 
tor function g(vvv • • • ) and integrates over velocity. The 
zeroth moment results from g = 1, the first from g = mv, 
the second from g = ^mvv, etc. The resulting equations 
are continuity equations for density, momentum, kinetic 
energy, etc. [52i l53l. l54| . 

Each momentum equation contains quantities that can 
be obtained from higher moments. Therefore the series of 
momentum equations has to be truncated in a convenient 
way, and one needs to introduce an assumption in the 
highest moment equation to be included. If we truncate 
after the first moment and simplify further [54[ in the 
momentum equation, we retain: 

^ + V-j e =5 (13) 
j e = -pEn e - DVn e (14) 

where n e is the electron density, j e = un e is the flux and 
u = (v) is the mean velocity of electrons. S is the source 
of electrons due to collisions and impact ionization, jj, 
represents the mobility and D is the diffusion matrix. 



The coefficients S, /i and D appearing in equations (fT3j) 
and (fl4"l) are to be obtained from elsewhere. One com- 
mon approach is to solve the Boltzmann equation (jT3J) 
for a homogeneous and constant electric field E within 
a background gas of constant density. In a uniform elec- 
tric field, the electrons gain energy from the field and 
loose it in inelastic collisions, reaching some steady state 
transport conditions [55l. |56| . By separating /(x, v, t) 
in i 7 '(v)n e (x, i), expanding F(v) in spherical harmon- 
ics, and truncating this series, a practically solvable set 
of ordinary differential equations results 57]. Solving 
these equation determines the electron velocity distribu- 
tion function f(E,v) from which [i(E) and D(E') with 
E = |E| can be calculated. Furthermore, the electron 
source term is written as 

S = K E| a(E), (15) 

when attachment and recombination can be neglected. 
Using these coefficients in a given gas and density as a 
function of the electric field is called the local field ap- 
proximation. 

Of course, this fluid model has to be extended by con- 
tinuity equations for other relevant excited or ionized 
species. For a nonattaching gas with neglected ion mobil- 
ity, the continuity equation for the density n p of positive 
ions has to be included 



Alternatively, the fluid model (fT3")) - (flT)]) can also be 
motivated b y p hysical considerations and conservation 
laws [II, [II II| . 

In the present paper, the highest possible consistency 
between particle and fluid model is achieved by determin- 
ing the transport coefficients and ionization rate fJ>(E), 
matrix D(E), and a(E) for the fluid model from swarm 
experiments in the particle model; this will be done in 
section III CI 

Solutions of the particle model and of the fluid model in 
local field approximation will differ when the electric field 
or the electron density vary rapidly in space or time as 
the electrons then will not fully "equilibrate" to the local 
electric field [22j. We will investigate these deviations 
further below. 



2. Numerical implementation 

The fluid equations are solved on a uniform grid where 
the electron densities n e and ion densities n p are calcu- 
lated at the centers of the grid cells. The densities can 
be viewed as averages over the cell like in the particle 
model. The field strength is also calculated at the cell 
centers. The electric field components are taken on the 
cell faces, where they determine the mass fluxes. 

The equations for the particle densities are discretized 
in space with the finite volume method, based on mass 
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balances for all cells. The particle densities are updated 
in time using the third order upwind-biased advection 
scheme combined with a two-stage Runge-Kutta method. 
For the details and the tests of the algorithm, we refer 

to a 

Analytical studies Qji, |36| and numerical investiga- 
tions H3] show that the ionization front is a pulled 
front; therefore a very fine numerical grid is required in 
the leading edge region of the front, and standard refine- 
ment techniques refining in the interior front region fail. 
Like for fronts in the particle model, we also have tested 
different numerical meshes for fronts in the fluid model, 
these fronts are treated in section IIII Al On a too coarse 
grid, the front moves too fast and is too smooth due to 
numerical diffusion of the electron density. To achieve the 
same numerical accuracy below 0.2 % as for the particle 
model, the fluid model requires an approximately four 
times finer mesh, namely dx — 0.575 /im and dt — 0.075 
ps. This mesh will be used below. 



C. Kinetics and transport of electron swarms in 
constant fields 

Swarm experiments deal with electron swarms mov- 
ing and multiplying in a constant electric field without 
changing it. If the field is high enough such that the 
electron number grows measurably, such a swarm is also 
called an avalanche. Swarms or avalanches in homoge- 
neous fields are an important experimental and theoret- 
ical tool to investigate the electron dynamics. 



1. Particle swarm kinetics: approach to steady state 

Particle swarm experiments can be used to determine 
transport coefficients, but also to study the particle ki- 
netics. We will first study the second issue, namely the 
relaxation of electrons to a steady state velocity distri- 
bution and the velocity distribution itself. Indeed the 
electron swarm will rapidly "equilibrate" to the applied 
field. In such a balanced state, the electrons on aver- 
age gain as much energy from the electric field as they 
loose in inelastic and ionizing collisions; this is how they 
reach an energy and velocity distribution specific for the 
electric field. The time that the electrons need to get 
in balance with the local electric field is an important 
indication for the validity of the local field approxima- 
tion. We therefore test it here within a particle swarm 
experiment. 

Fig. [3] shows how different electron swarms converge 
to the same mean energy within a field of —100 kV/cm. 
The experiment starts with a group of electrons of identi- 
cal velocity directed in the electron drift direction; their 
kinetic energy is 50, 5 and 0.5 eV, respectively. When 
the swarms start to drift, their average energies converge 
to the same constant value within at most 3 ps in all 
three cases. Here the average energy is used as the sim- 
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Figure 3: Average electron energy as a function of time for 
three different electron swarms in a field of —100 kV/cm. 
Starting with electron swarms moving into the drift direction 
with an identical kinetic energy of 50 eV (solid), 5 eV (dotted) 
and 0.5 eV (dashed), all electron swarms approach a mean 
electron energy characteristic for the applied field within 2 to 
3 ps. 



plest indication of their energy and velocity distribution 
function. 

After reaching steady state motion, swarms demon- 
strate a steady state distribution of electron velocities 
and energies in a given electric field. Fig. 2] shows the 
distribution of the longitudinal electron velocity v z in 
an electric field of 50, 100, 150, and 200 kV/cm. The 
figure shows that with increasing field, the electron ve- 
locity distribution deviates more and more from the 
Maxwellian profile and therefore from symmetry about 
velocity v z = 0; rather an increasing number of electrons 
flies in the direction of the field and a decreasing umber 
against it, and the number of electrons with high kinetic 
energy increases. 



2. Gaussian swarm profiles and transport coefficients 

Swarm experiments are used as well to determine mo- 
bilities, reaction rates and diffusion constants experimen- 
tally [59( | . An electron swarm drifts, broadens and grows 
under the influence of a constant electric field. The same 
experiment is performed here for this purpose, but now 
numerically with the particle model. 

We here recall the essentials: A single electron will 
generically evolve into a swarm that has a Gaussian pro- 
file in space. In terms of the fluid model (fT3"]) - (fTuT) . this 
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Figure 4: Distribution of the electron velocity v z in the lon- 
gitudinal direction in particle swarm experiments in fields of 
50, 100, 150, and 200 kV/cm. The higher the field, the more 
the distribution deviates from Maxwellian and from symme- 
try about velocity v z — 0. 



Gaussian distribution is given by 

-(x 2 + y 2 )/(4irD T (t-t Q )) 



n e (x,y,z,t) oc e 



fj,Ea(E) t e _ 



4nD T (t-t ) 
\z - z - nEt) 2 /{4*D L (t - t Q )) 



^4TTD L (t-t ) 



; (17) 



(The dischar ge s pecific context of this solution can be 
found in [58L l60j|.) Here the center of the package is at 
(x, y, z) = (0, 0, zq) at time t — 0, and the field is in the z- 
direction. The longitudinal and transversal components 
of the diffusion matrix are denoted as Dl and Dt- 

The transport and reaction coefficients can be deter- 
mined from this equation by 



KE)\E\ 
a(E) 
D T {E) 
Dl(E) 



*2 - *1 

1 lnJV e (t 2 )-lniVe(ti) 



H(E)\E\ t 2 -h 
{x\t 2 )+y\t 2 ))-{x^t 1 ) + y^t l )) 



4(t a -*i) 

{[z{t 2 )-{z{t 2 ))Y)-{[z{t 1 )-{z{t 1 ))Y) 



(18) 
(19) 
(20) 



2(t 2 - tl) 



(21) 



where the spatial averages of observables O are taken 
over the electron density (flT|) 



(0(t)) = J dx dy dz O n e (x, y, z, t) for the fluid model, 

(22) 



and N e (t) is the total number of electrons at time t. 

When swarm experiments are performed in the particle 
model, transport and reaction coefficients are determined 
by (fT8|) -([2T |) as well, and averages are now performed 
over all electrons 

(£?(*)) = ——t °i( t ) for thc P article model, 

3 — 1 

(23) 

in agreement with equation ([3]). This is how the coeffi- 
cients for the fluid model will now be determined from 
the particle model. 



3. Fluid parameters determined from particle swarms 

We have determined n(E), D T (E), D L (E), and a(E), 
and also the average electron energy e(E) in particle 
swarm experiments for 42 different background electric 
fields ranging from 2 kV/cm to 205 kV/cm. 

To obtain the transport coefficients and mean values 
with satisfactory statistics, one needs a sufficient number 
of electrons that have experienced an adequate number 
of collisions. The experiments start from a number of 
electrons at the same position (i.e. located in a single 
point, which is, a Gaussian with zero width), and end 
with a swarm of electrons with a Gaussian distribution 
as described in equation (fTT)) . Because the ionization 
rate depends strongly on the electric field strength, the 
number of initial electrons and the simulation time is cho- 
sen according to the fields. For example, the simulation 
starts with 10 6 electrons at 2 kV/cm and lasts for 1500 
ps, but for 205 kV/cm, the simulation starts with 10 2 
electrons and ends with 4 • 10 6 electrons after 30 ps. 

As there is some initial transient during which thc elec- 
trons "equilibrate" to the field and approach a Gaussian 
density profile, the transport and reaction coefficients are 
evaluated according to equations (fT8|) - (f2~Tj) at appropri- 
ate times ti 2. We choose t 2 as the end of a swarm exper- 
iment, and t\ = t 2 /2 in the middle of an experiment. In 
view of the relaxation times below 3 ps evaluated above, 
this choice of t\ is on the very safe side. 

The numerical results for /i, Dt, Dl, a and e for differ- 
ent electric fields E = |E| are presented in Fig.[5]together 
with empirical fit formulas. These formulas are 
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Figure 5: Electron mobility, diffusion rates, ionization rate and average electron energy in nitrogen. Plotted are the reduced 
coefficients fj,p, a/p, Dt, Dl and e as a function of reduced field E/p at room temperature. Our units are related to other 
commonly used units like I kV cm -1 bar -1 = 1.316 V cm -1 torr -1 = 0.424 Townsend at T=300 K. 
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e(E) = eV • exp 
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, (28) 



under normal conditions. 

These coefficients will be used in the fluid model (|13|) 
(|16p to reach optimal agreement between particle and 
fluid model. 



III. SIMULATIONS OF PLANAR FRONTS 

A. Concepts and set-up of planar ionization fronts 

1. The role of planar fronts in the inner analysis of the 
streamer structure 

Fluid model simulations of streamers within the past 
20 years (see e.g., [U El, El HI) have shown that the 
streamer head is surrounded and preceded by an ion- 
ization front that propagates into the non-ionized gas. 
Within the ionization front, ionization grows until space 
charge effects set in. The formed space charge layer is 
much thinner than the radius of the streamer, it leads 
to a screening of the electric field in the interior of the 
streamer head and to a field enhancement ahead of it. 
Therefore the field dependent ionization reaction coeffi- 
cient a(E) is enhanced ahead of the space charge layer 
and suppressed behind it. The space charge layer around 
the streamer is shown in Fig. [I] f° r a further discussion of 
the three-dimensional structure and growth of streamers, 
we refer to the literature, see, e.g., [27] . 

It is clear that the full configuration of the electric field 
can only be analyzed within a two- or three-dimensional 
setting. On the other hand, within the inner structure of 
ionization front and space charge layer, the electric field 
does not deviate much from a planar configuration. To 
analyze the processes within the ionization front in de- 
tail, it is therefore advisable to study the inner structure 
of a planar front. This will be done here. The results 
can be put in further physical context through a sep- 
arate analysis on the inner and the outer scale of the 
structure as commonly done in hydrodynamic boun dary 
layer analysis, reaction-diffusion systems etc. (39l.l40l[4l|. 



by imposing periodic boundary conditions at the lateral 
boundaries. Furthermore, the electric field is calculated 
only in the forward direction z through 



E z (z,t) = E z (z ,t)+ dz' 

J z J A 



dx dy e(n. p — n e )(x, y, z', t) 



A 



'29; 

which is the one-dimensional version of ([2]). This means 
that fluctuations of the transversal field due to density 
fluctuations in the transversal direction are not included. 
In fact, the numerical implementation of (f2"9"|) is per- 
formed on a grid in the forward direction only as dis- 
cussed in section Hi A 21 

The density fluctuations projected onto the forward 
direction depend on the transversal area A over which the 
averages are taken. When A increases, the total number 
of electrons in the simulation increases proportionally to 
A, while the relative density fluctuations decrease like 
\j\f~A. Therefore some intermediate value of the area 
A has to be chosen: On the one hand, there should be 
a sufficient number of electrons to reach a satisfactory 
statistics, but on the other hand, there shouldn't be so 
many electrons that the computer runs out of memory 
within the time interval of interest. In the simulation, 
we use a small A for high electric field and a large A fol- 
low electric field. For example, we choose the transversal 
averaging area as A = 6 At x 6 At for —100 kV/cm, 
but for —50 kV/cm, the transversal averaging area is 
A = 20 A£ x 20 Ai, here Ai = 2.3 fim. 

In the z direction, the system length is 500 At which 
allows the front to propagate freely for all runs reported 
in this paper. The electric field in the non-ionized region 
at large z is specified by E — E + z, where E + < and z 
is the unit vector in the z direction. 

In the simulations, two different initial conditions are 
used. In the first, N e electrons are evenly distributed in 
a thin layer of area A with an extension of 19.5 At < z < 
20.5 At in the field direction. Choosing 



(30) 



N e /A= \E+\ •e /e; 



the field behind the layer is screened according to 
Another choice is to begin with a few seed electrons which 
will create an ionization avalanche and form a charge 
layer later. 



B. Planar fronts in the particle model 



2. Construction of planar fronts in the particle model 



Qualitative discussion of typical results 



The construction of a planar front is straight forward 
in the fluid model: gradients V are simply evaluated in 
one spatial direction. We choose this direction to be the 
z direction. In the particle model, on the other hand, 
electrons move in all three spatial dimensions. There- 
fore a three-dimensional setting has to be retained. An 
essentially one-dimensional setting is achieved by consid- 
ering only a small transversal area A of the front and 



We first present results in a field of E + = —100 kV/cm. 
The initial condition is a thin electron layer with total 
electron number N e as in (f3"0"|) , screening the electric field 
behind the layer. Fig. [5] presents the evolution at times 
t\ = 450 ps (left) and £2 = 900 ps (right). Panel (a) shows 
the density distribution of the electrons (solid) and the 
ions (dashed). Panel (b) shows the net negative charge 
distribution (solid) and the electric field (dashed). Panel 
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Figure 6: Spatial profile of a particle front within a field of —100 kV/cm at time ti = 450 ps (left) and t% = 900 ps (right), 
(a) Electron density (solid) and ion density (dashed) , (b) electric field (dashed) and net charge (solid) , (c) density of electrons 
with energy above eV, 15.6 eV, 20 eV, 30 eV and 50 eV, (d) zoom into the front region with high electron energies, same 
quantities as in (c). 



(c) shows the total charge density of electrons and the 
charge density of electrons with an energy higher than 0, 
15.6, 20, 30, and 50 eV, where 15.6 eV is the ionization 
energy. Panel (d) zooms into panel (c), both in space 
and in charge densities. 

The figure shows an ionization front propagating to 
the right; up to fluctuations, the spatial profiles are es- 
sentially unchanged, therefore the front velocity v is es- 
sentially constant as well. The front carries an overshoot 
of electrons, generating a thin space charge layer that 
screens the electric field behind the front. In this screened 
interior region, the electron and ion density reach an 
equal constant density n~ = n~ . (Upper indices "±" 
indicate quantities before "+" and behind "— " the ion- 
ization front.) The qualitative features of the front are 
the same as those in the fluid model [HI]. 

Electrons with energies above the ionization threshold 
of 15.6 eV are shown in the lower two panels in Fig. [51 



they exist essentially only in the high field region. Elec- 
trons with energies above 30 eV are so rare that they are 
hardly seen even on the scale of panel (d) . Electrons with 
energy above 50 eV exist, but cannot be distinguished 
within this plot. The profiles of high energy electrons 
also move with the whole front without change of shape, 
up to fluctuations. 



Following the track of single electrons of high energy, 
we found that they gain and loose energy in few collisions 
within a few ps and do not run away. In that sense 
they are in a fast dynamic equilibrium with the electrons 
at lower energies. This observation agrees with the fast 
relaxation of 50 eV electrons travelling in the forward 
direction whose energy relaxation is shown in figure [31 
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Table I: Numerical results on planar fronts in the particle Table II: Exact analytical results for planar fronts in the fluid 
model: front velocity v and ionization level n~ behind the model, evaluated with the transport coefficients from ((24]) — 
front as a function of the electric field E + ahead of the front. (??): front velocity v* and electron density decay length I* 

as a function of the electric field E + . 



2. Quantitative results in different fields 

The ionization front in a given field E + is characterized 
by a velocity v, a degree of ionization n~ — rip behind 
the front and an electron energy distribution in the high 
field region. We now present these quantities in detail. 

We define the front position as the position of the max- 
imal electron density. Table [J summarizes our numerical 
results on the front velocity v as well as on the saturated 
electron density n~ behind the front as a function of the 
electric field E + immediately ahead of the front. 

C. Planar fronts in the fluid model 
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Table III: The degree of ionization n^T behind the front in the 
fluid model as a function of field E + : numerical result n^T and 
analytical upper bound n~ bound as derived in the appendix. 



For planar fronts in the fluid model, there are not only 
numerical, but also analytical results, both agree within 
the numerical accuracy. First, the velocity of the front 
in a field E + is given by 



v* = (x(E+) \E+\ ± 2^D L (E+) f i(E+) \E+\ a(E+), 

_ (31) 
according to [IS, [38| with a slight generalization along 
the general lines discussed in [3g|. Note that for the ini- 
tial conditions treated in this paper, the velocity v* is 
approached from below after an algebraically slow relax- 
ation 36]. 

The electron and ion densities decay exponentially like 
e -z/e j n fo e leading edge of the front where the electric 
field is approximately E + [TH, [3(| ■ The decay length is 



D L (E* 



H(E+) \E+\ a(E+)' 



(32) 



These analytical results are summarized in Table [TTJ 

For the degree of ionization behind the front n~ , there 
is no closed analytical solution. The ionization behind 
the front can be derived numerically. Furthermore, in 
the appendix, we derive an analytical upper bound for 
this quantity, namely 



< 



e 



E+\ 



a(e) de. 



(33) 



Numerical result and analytical bound are summarized 
in table 1TTT1 



D. Comparison of planar fronts in particle and 
fluid model 



1. Detailed investigation in a field of — 100 kV/cm 

Now the stage is set to compare planar fronts in the 
particle model to those in the fluid model, The compari- 
son is first done in detail for a planar front propagating 
into a field of —100 kV/cm. Both fluid and particle sim- 
ulations are carried out in the same setup starting from 
the same initial conditions, i.e., from an electrically neu- 
tral group of 100 electrons and ions evenly distributed 
within the thin layer 19.5 A£ < z < 20.5 A£ and within 
the transversal area A = 6 A£ x 6 A£. 

Fig.[7]shows the temporal evolution of the planar front. 
We compare the spatial profile of the electron density 
(solid line) and ion density (dotted) in a particle simu- 
lation with the electron density (dashed) and ion den- 
sity (dot-dashed) in a fluid simulation. Two features are 
clearly visible: First, the particle and the fluid front move 
with approximately the same velocity and the same den- 
sity profile in the leading edge of the front where the 
electric field does not vary yet. Second, the maximal 
electron density in the front and the saturation level of 
the ionization behind the front in the particle model are 
about 20 % higher than in the fluid model. These results 
of visual inspection agree with those of tables I and III 
for a field of -100 kV/cm. 

As we have excluded other reasons of the discrepancy 
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Figure 7: Temporal evolution of the electron and ion densities 
in a planar front in a field of 100 kV/cm. Shown are the spatial 
profiles of electron and ion densities derived with the particle 
or the fluid model at time steps t = ns, 0.09 ns, 0.18 ns, 
0.27 ns, . . ., 0.9 ns (solid lines: n e , par t, dashed: n e ,fiuid, dot- 
ted: n PtPa rt, dot-dashed: n v ji ui d). 
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Figure 8: Zoom into the particle ionization front of Fig. [7]in a 
field of —100 kV/cm at time 0.63 ns: shown are electron den- 
sity distribution n e (solid line), local average electron energy 
£ P art (dotted line), local average electron energy elfa accord- 
ing to the local field approximation (dashed line) , and electric 
field strength E (dot-dashed line). 



like numerical discretization errors or inconsistent trans- 
port and reaction coefficients, deviations must be due 
to the approximations in the fluid model, and a closer 
inspection shows that we should focus on the electron 
energies. 

Fig. [8] zooms into the ionization front shown in Fig. [7] 
at time t = 0.63 ns. Here we show the electron density 
(solid line) and electric field (dot-dashed line) in the par- 
ticle model. Furthermore the local mean energy of the 
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Figure 9: The electron energy distribution function is mea- 
sured as P(e)/y/e, where P{s) is the probability of electron 
energy e. Here we show the electron energy distribution at the 
front region of Fig. |8]where E w —90 kV/cm (solid line) and 
in the swarm experiments in constant fields of —90 kV/cm 
(dashed line) and —100 kV/cm (dotted line). 



electrons in the particle model is indicated with a dot- 
ted line. Finally, the mean electron energy according to 
the local field approximation e(E) is derived from the 
local field E and equation (??) for e(E); it is indicated 
with a dashed line. It can be seen that the average elec- 
tron energy nicely follows the local field approximation in 
the interior of the ionized region while it is considerably 
higher in the region where the electric field is large and 
the electron density decreases rapidly. 

This deviation is analyzed in more detail in Fig. [5] 
where not only the mean energy, but the full electron 
energy distribution is shown. This is done for the par- 
ticular spatial region of the front where the electric field 
has decreased by 10 % to a value of —90 kV/cm. More 
precisely, electrons are collected from the first cell where 
\E\ > 90 kV/cm, searching with increasing z. The av- 
erage field in these cells is about E = —91.47 kV/cm. 
To reach a satisfactory statistics, the electrons are col- 
lected from 10 different instants of time with a temporal 
distance of 30 ps between the two consecutive sampling 
times to ensure statistical independence. While these 
data are plotted as a solid line, for comparison the elec- 
tron energy distributions in the swarm experiments in a 
constant field of —90 and —100 kV/cm are plotted as a 
dashed or dotted line, respectively. Analyzing the elec- 
tron energy distribution at low energies e < 10 eV, the 
ionization front at a local field around — 90kV/cm and 
the swarm experiment in a constant field of —90 kV/cm 
are quite similar while the distribution for a swarm in a 
field of —100 kV/cm is clearly lower. On the other hand, 
for electron energies above 20 eV, the energy distribution 
in the ionization front at a local field of —90 kV/cm actu- 
ally lies closer to the distribution of the swarm at —100 
kV/cm than to that at —90 kV/cm. This observation 
not only confirms that the average electron energy in the 
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Figure 10: Relative difference of ionization level behind the 
front (|34p (dotted, above), front velocity (|35p (dashed, be- 
low), and mean electron energy in the leading edge (|36[) (solid, 
middle) between particle and fluid model as a function of the 
electric field E + ahead of the front. 



leading edge of the ionization front is higher than in the 
local field approximation, as is consistent with Fig. [HI but 
it indicates that the higher mean energies correspond to 
a higher population of electron energy states above the 
ionization threshold. We therefore expect that also the 
ionization rates are higher in the particle model than in 
the local field approximation. 



2. Results for other fields 

Having analyzed the front propagating into a field of 
E + = —100 kV/cm in detail, we now summarize equiv- 
alent results for fields ranging from —50 to —200 kV/cm 
in Fig. [TU] The figure shows the relative difference 



'e, fluid 



(34) 



'e, fluid 



of the saturated electron density behind the front be- 
tween the particle and the fluid model, the relative dif- 
ference 



Wpart — Wfluid 
Wfluid 

of the front velocity and the relative difference 

(gpart) ~ (£LFA) 

(elfa) 



(35) 



(36) 



of the mean electron energy between the particle model 
and the local field approximation at a point in the ioniza- 
tion front where electron densities are low and the electric 
field is just slightly screened to the value E = 0.98 E + . 
All differences increase with increasing field. 



E. Interpretation of results 

1. Electric field gradient, electron motion in the front and 
electron energy overshoot 

A few authors have dealt with particle discharge mod- 
els with strong gradients both in the electric field and in 
the particle density. They analyzed deviations between 
particle models and the fluid model in local approxima- 
tion by means of the Boltzmann equation, as treated in 
section III B 11 They analyzed the case of stationary gra- 
dients of the electric field (5 5| . [5 6Tl and the case of positive 
streamer ionization fronts|22| . They suggested correc- 
tions to the fluid model both due to field gradients and 
to density gradients. 

The situation in negative streamer ionization fronts is 
somewhat different: a negative ionization front in nitro- 
gen moves approximately with the electron drift velocity 
determined by the electric field in the leading edge of the 
ionization front. This means that the electrons in this 
leading edge region on average move within a stationary 
electric field. More precisely, the front velocity v*{E + ) 
from (f3"Tj) is slightly larger than the electron drift velocity 
n{E + ) \E + 1 in the electric field ahead of the front. This is 
because the front is not carried by electron drift only, but 
also by diffusion and creation of new free electrons. In a 
coordinate system moving with the ionization front, the 
existing electrons therefore on average move backwards. 
This motion is the faster, the further they are behind. 

Fig. [5] shows that the high energy tail of the electron 
energy distribution in the ionization front at the position 
where the field is —91.47 kV/cm, is closer to the swarm 
experiment at —100 kV/cm than to that at —90 kV/cm. 
One could interpret these results by assuming that the 
electrons have gained their energy distribution in a field 
of close to —100 kV/cm and then are transported back- 
wards to where the field is only —90 kV/cm. 

However, Fig. [SJ shows that this interpretation cannot 
be true. If the electron energy distribution would first 
"equilibrate" to a field of —100 kV/cm and then par- 
tially be transported backwards, then the mean electron 
energy (e pa rt) would everywhere be below the mean elec- 
tron energy (elfa) at —100 kV/cm in local field approxi- 
mation. But clearly there is an electron energy overshoot 
in the leading edge of the ionization front. So the most 
prominent deviation from the fluid model visible in Fig. [5] 
occurs in the region were the electric field is almost con- 
stant; furthermore, if it would be a consequence of the 
local field approximation, the deviation would have the 
opposite sign. 



2. Density gradient and relation between front and swarm 
experiments 

We conclude that the electron energy overshoot in the 
leading edge of the ionization front has to be due to the 
electron density gradient rather than to the field gradient. 
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The effect of density gradients can be tested in swarm 
experiments (cf. section lTl CP in a constant electric field as 
well. It turns out that such a test goes beyond qualitative 
results and actually allows for a quantitative comparison 
of fronts and swarms as will be explained below. 

The key to the quantitative comparison is based on 
two closely related facts: (i) the density profiles in the 
leading edge of a swarm or avalanche in a constant field 
and in the leading edge of an ionization front penetrating 
the same field are both given by 

n e ~e- z / r , (37) 

where t* is given by ([3^)) and in table II. (ii) Also the 
velocities of swarm and front have the same value v* from 
(|3T|) and table II. 

This relation between velocities and decay rates of 
swarms and fronts holds generally for any so-called pulled 
front whose velocity is determined in the linearly unsta- 
ble region ahead of the front, as is discussed in a general 
setting in j2f|, specifically in section 2.5.1. The state- 
ment can be verified on the explicit form of the Gaussian 
profile (jTTJ) : in the coordinate 



Zq — v*t 



(38) 



moving with velocity v* , the Gaussian swarm distribution 
(fT7|) can be written as 



, fi\E\at 



-{z-z -^Et) 2 /(AD L t) 

s/4nD L t 



,-e/(iD L t) 



y/4irD L t 
(39) 

We now test the above theoretical predictions on the 
particle model for a swarm and a planar front in a field 
of —100 kV/cm. The results are shown in Fig. [TTJ It 
should be remarked that our electron swarm has a Gaus- 
sian density distribution both in the longitudinal and 
in the transversal direction. To focus on the profile in 
the longitudinal direction, the density in the swarm in 
Fig. [TTJa is taken as the density on the longitudinal axis. 
Fig. [TTJa shows the electron density profile of the swarm 
decaying at both edges and the profile of the planar front 
that grows and saturates to a constant level. As the den- 
sities are plotted on a logarithmic scale, an exponential 
decay like in (|3T[) amounts to a straight line with slope 
— l/£* in the plot. Despite density fluctuations and slow 
transients in the buildup of the profile [36[, Fig. [TTJa in- 
deed shows that swarm and front have a similar decay 
profile on the right hand side. 

This sets the stage to compare now the electron en- 
ergies of swarm and front in Fig [TTJb. The dotted and 
the dashed line reproduce lines from Fig. [5] for the pla- 
nar particle front, they show the actual mean electron 
energy (e pa rt) (dotted) and the mean electron energy ac- 
cording to the local field approximation (elfa) (dashed). 
For the electron swarm, the mean electron energy along 
the swarm axis is indicated as a solid line in Fig [TTJb. 
It is remarkable that this energy within the swarm is 
not constant in space though the electric field has a con- 
stant value. Rather the energy increases almost linearly 



along the axis. While the average energy of the swarm 
is around 8.5 eV, producing the respective value for the 
local field approximation e(E) in Fig. 5(e)[ the average 
energy at the leading edge of the swarm reaches values 
above 10 eV. 

Now the leading edge regions of swarm and front need 
a closer inspection. The leading edge is defined as the 
region where the electric field is (almost) constant and 
where the electron densities both approach the profile 
(|37|) . It is in this region where the electrons attain the 
highest mean energies, and FigQTJb shows that the mean 
energy profiles of swarm and front in this region are vir- 
tually identical up to fluctuations. We conclude that the 
high electron energies within the leading edge of the ion- 
ization front that are shown in figures [5] and [5J are due 
to the electron density gradient, and can be studied in 
the leading edge of a swarm experiment as well. 



3. The ionization density behind the front 

The discussion above shows that within the electron 
density gradient, the fast electrons are on average a bit 
ahead of the slower electrons as they have a higher av- 
erage velocity in the forward direction. This is under- 
standable since on the one hand, each single fast electron 
looses its energy over a length of ~ 1 fj,m (which is the 
electron drift velocity fi(E)E times the relaxation time 
~ 2 ps from figure^]), but on the other hand, the gradient 
length I* « 2.7 ^m of the density profile is of the same 
order. Clearly, yhis effect becomes the more pronounced, 
the higher the field. 

Now the electron energies in the leading edge substan- 
tially exceed the local field approximation and indicate 
the presence of a larger fraction of electrons with energies 
above the ionization threshold. This leads also to higher 
ionization rates than estimated by the local field approx- 
imation ol[E) (??). As the ionization level behind the 
front n~ = n~ is well approximated by the effective ion- 
ization rate a integrated over the fields within the front 
([3"3"| (recall also the argument in the appendix), it is clear 
that the ionization level behind the front is higher in the 
particle model than in the fluid model. 



IV. CONCLUSION 

Negative streamer ionization fronts in pure nitrogen 
were investigated in planar approximation, both follow- 
ing the kinetics of single electrons in a particle model, 
and in a fluid model in local field approximation. As pa- 
rameter functions for the fluid model were derived from 
swarm experiments in the particle model and numerical 
errors are under strong control, a discrepancy between 
results of particle and fluid model must be attributed to 
the approximations made in the fluid model. 

For electric fields immediately ahead of the front of 50 
kV/cm or lower — the statements hold for nitrogen un- 
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Figure 11: Comparison of an electron swarm and a planar ionization front, both in the particle model and in a field of —100 
kV/cm. (a) Electron density on a logarithmic scale for swarm (solid) and front (dotted) as a function of z. (b) Mean electron 
energy in swarm (solid) and front (dotted) and in the front assuming the local field approximation (dashed) as a function of z. 



der normal condition, they are easily extended to other 
densities by introducing the reduced electric field E/N — 
, particle and fluid model essentially agree in front speed, 
profile and ionization level behind the front. When the 
field increases, the velocity does not vary much between 
both models, but the ionization level behind the front 
is substantially larger, as is demonstrated for a field of 
100 kV/cm in figure [7] In figure ITD1 the differences be- 
tween particle and fluid model are summarized for the 
full field region investigated. 

Can the discrepancy between the models be attributed 
to a particular physical mechanism and to a particular 
part of the front region? In figure [8] we find the largest 
discrepancy between fluid and particle model in the lead- 
ing edge of the front where the electron density is very 
low and where the electric field is hardly screened. Here 
the electrons have an average energy considerably higher 
than estimated by the local electric field. This effect can 
only be explained by an effect of the electron density 
gradient, not of the electric field gradient. Indeed, an 
electron swarm or avalanche in a constant electric field 
has the same velocity and the same density gradients in 
its leading edge and shows the same electron energy over- 
shoot in its front part, as is illustrated in figure [TT1 In 
essence, in high fields the gradient length of front and 
avalanche £* becomes of the order of the electron energy 
relaxation length, therefore the fast and energetic elec- 
trons can get ahead of the slower ones. 



These higher electron energies in the leading edge of 
the front lead to higher ionization rates than the ion- 
ization rate a(E) in local field approximation — indeed, 
a(E) is the total effective ionization rate of a complete 
electron swarm (as derived in section lTl CI) while the elec- 
trons in the leading edge of the swarm have higher ion- 
ization rates. 

Now the ionization level behind the front can be ap- 
proximated by integrating the effective ionization rate 
a(E) over the electric field strength E, independently of 
the precise spatial structure of the front. This result is 
derived in the appendix. This implies that higher elec- 
tron energies and higher effective ionization rates in the 
leading edge of the front directly lead to higher ioniza- 
tion levels behind the front, even though only very few 
electrons are involved in this dynamics in the low density 
region. 

The effect of few electrons of high energy is much more 
severe for the ionization rates a than for the average drift 
velocity iiE, therefore the most pronounced effect is seen 
in the ionization levels behind the front, and less in the 
front velocities. 

We finally remark that next to explicit predictions, our 
work has delivered two useful insights: (i) the physical 
discrepancies between particle and fluid model lie in the 
leading edge of the front, though the effect is not so much 
seen in the velocity, but much more in the ionization 
level behind the front. This gives a clue for a numerical 
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strategy combining efficient features of fluid and particle 
model, (m) The essential features in the leading edge of 
the front are equally present in the leading edge of an 
electron swarm or avalanche in a constant field, where it 
can be studied much easier. 



while keeping the equations unchanged otherwise. This 
change of the source term means that more ionization 
is created at each part of the front, and therefore the 
ionization level behind the front will be higher. This is 
because in a negative front, drift current and diffusion 
current point into the same direction, and therefore 
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-U{E)-Vn e -/j(E) En e > fj,(E) |E| n e . (A3) 

The dominant part of the current is the drift or Ohmic 
part, therefore the results with the changed source term 
(|A2|) should be close to the original problem. 

According to Maxwell's equations, the divergence of 
the total current vanishes 



Appendix A: APPROXIMATING THE 
IONIZATION LEVEL BEHIND THE FRONT 

In previous papers [H|,[6lj], we have established that for 
a planar, uniformly propagating front with field indepen- 
dent electron mobility fi(E) = const, and for vanishing 
electron diffusion D(E) = 0, there is a unique relation 
between the field ahead of the front E + , the ionization 
rate function a(E) and the degree of ionization n~ = n~ 
behind the front, namely 

eo r lE+l 

n~ = n~ = — / a(e) de (Al) 
e Jo 

in dimensional units. 

Such a simple identity does not hold for the full fluid 
model @, (|13[) ~ (fTr?)) with electron diffusion, but we can 
establish the expression (|A1[) as an upper bound for the 
free electron density n~ behind the front; this upper 
bound is actually a very good approximation to the real 
value as tablc lllll shows. The upper bound is constructed 
as follows: 

In the ionization source term S (fT5} . we follow the 
usual procedure to take only the drift term of the elec- 
tron current into account. We note at this place that 
this assumption requires some reconsideration, and one 
could argue that the complete current should be taken 
into account and that a decomposition into a drift and a 
diffusion part is artificial. This argument will have to be 
elaborated at some other place. We here just change the 
source term into 

Si = Ue| <*(E), (A2) 



V-(e 3 f E+j)=0, (A4) 

where j is the electric current; in our case it is j = — e j e 
with j e from (fl4|) . If the front is planar, and if the electric 
field E + in the non-ionized region does not change in 
time, then it follows immediately that 

e d t E = e j e (A5) 

Now the growth of the ion density (|13|) is bounded by the 
source term (| A2|) and the identity (|A4[) is inserted: 

d t n p = S = n(E) |E| n e a{E) 

<S j =j e a(E) = ^d t Ea(E). (A6) 

The first and last expression in this inequality can be 
integrated in time with the result 

e /' B(x ' <) 

n p (x, t) — n p (x, 0) S — / a(e) de. (A7) 

~ e Je(x,0) 

Now the time interval [0, t] is taken in such a way that 
it contains the time range in which the ionization front 
passes over the point x of observation. n p (x,0) is then 
the ion density before and n p (x, t) the ion density behind 
the front. As a result, we find that the expression (|A1[) 
indeed is an upper bound and good approximation for 
the ionization level behind the front. 

We finally note if there is electron attachment and pos- 
itive and negative ions n Pi „ are formed, the statement 
stays true for the total ion charge density n = n p — n n if 
the total source and sink term for the ion density n can 
be written in the form (|A2[) . 
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